home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / ode-initval / gear2.c < prev    next >
Encoding:
C/C++ Source or Header  |  2002-04-18  |  5.6 KB  |  217 lines
  1. /* ode-initval/gear2.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. /* Gear 2 */
  21.  
  22. /* Author:  G. Jungman
  23.  */
  24. #include <config.h>
  25. #include <stdlib.h>
  26. #include <string.h>
  27. #include <gsl/gsl_math.h>
  28. #include <gsl/gsl_errno.h>
  29. #include "odeiv_util.h"
  30. #include <gsl/gsl_odeiv.h>
  31.  
  32.  
  33. /* gear2 state object */
  34. typedef struct {
  35.   int primed;                /* flag indicating that yim1 is ready */
  36.   double t_primed;           /* system was primed for this value of t */
  37.   double last_h;             /* last step size */
  38.   gsl_odeiv_step * primer;   /* stepper to use for priming */
  39.   double * yim1;             /* y_{i-1}    */
  40.   double * k;                /* work space */
  41.   double * y0;               /* work space */
  42.   int stutter;
  43. }
  44. gear2_state_t;
  45.  
  46. static void *
  47. gear2_alloc (size_t dim)
  48. {
  49.   gear2_state_t *state = (gear2_state_t *) malloc (sizeof (gear2_state_t));
  50.  
  51.   if (state == 0)
  52.     {
  53.       GSL_ERROR_NULL ("failed to allocate space for gear2_state", GSL_ENOMEM);
  54.     }
  55.  
  56.   state->yim1 = (double *) malloc (dim * sizeof (double));
  57.  
  58.   if (state->yim1 == 0)
  59.     {
  60.       free (state);
  61.       GSL_ERROR_NULL ("failed to allocate space for yim1", GSL_ENOMEM);
  62.     }
  63.  
  64.   state->k = (double *) malloc (dim * sizeof (double));
  65.  
  66.   if (state->k == 0)
  67.     {
  68.       free (state->yim1);
  69.       free (state);
  70.       GSL_ERROR_NULL ("failed to allocate space for k", GSL_ENOMEM);
  71.     }
  72.  
  73.   state->y0 = (double *) malloc (dim * sizeof (double));
  74.  
  75.   if (state->y0 == 0)
  76.     {
  77.       free (state->k);
  78.       free (state->yim1);
  79.       free (state);
  80.       GSL_ERROR_NULL ("failed to allocate space for y0", GSL_ENOMEM);
  81.     }
  82.  
  83.   state->primed = 0;
  84.   state->primer = gsl_odeiv_step_alloc (gsl_odeiv_step_rk4imp, dim);
  85.   state->last_h = 0.0;
  86.  
  87.   return state;
  88. }
  89.  
  90. static int
  91. gear2_apply(void * vstate,
  92.             size_t dim,
  93.             double t,
  94.             double h,
  95.             double y[],
  96.             double yerr[],
  97.             const double dydt_in[],
  98.             double dydt_out[],
  99.             const gsl_odeiv_system * sys)
  100. {
  101.   gear2_state_t *state = (gear2_state_t *) vstate;
  102.  
  103.   state->stutter = 0;
  104.  
  105.   if(state->primed == 0 || t == state->t_primed || h != state->last_h) {
  106.  
  107.     /* Execute a single-step method to prime the process.  Note that
  108.      * we do this if the step size changes, so frequent step size
  109.      * changes will cause the method to stutter. 
  110.      * 
  111.      * Note that we reuse this method if the time has not changed,
  112.      * which can occur when the adaptive driver is attempting to find
  113.      * an appropriate step-size on its first iteration */
  114.  
  115.     int status;
  116.     DBL_MEMCPY(state->yim1, y, dim);
  117.     
  118.     status = gsl_odeiv_step_apply(state->primer, t, h, y, yerr, dydt_in, dydt_out, sys);
  119.  
  120.     /* Make note of step size and indicate readiness for a Gear step. */
  121.  
  122.     state->primed = 1;
  123.     state->t_primed = t;
  124.     state->last_h = h;
  125.     state->stutter = 1;
  126.  
  127.     return status;
  128.   } else {
  129.     /* We have a previous y value in the buffer, and the step
  130.      * sizes match, so we go ahead with the Gear step.
  131.      */
  132.  
  133.     double * const k  = state->k;
  134.     double * const y0 = state->y0;
  135.     double * const yim1 = state->yim1;
  136.  
  137.     const int iter_steps = 3;
  138.     int status = 0;
  139.     int nu;
  140.     size_t i;
  141.  
  142.     DBL_MEMCPY(y0, y, dim);
  143.  
  144.     /* iterative solution */
  145.  
  146.     if(dydt_out != NULL) {
  147.       DBL_MEMCPY(k, dydt_out, dim);
  148.     }
  149.  
  150.     for(nu=0; nu<iter_steps; nu++) {
  151.       int s = GSL_ODEIV_FN_EVAL(sys, t + h, y, k);
  152.       GSL_STATUS_UPDATE(&status, s);
  153.       for(i=0; i<dim; i++) {
  154.         y[i] = ((4.0 * y0[i] - yim1[i]) + 2.0 * h * k[i]) / 3.0;
  155.       }
  156.     }
  157.  
  158.     /* Estimate error and update the state buffer. */
  159.     for(i=0; i<dim; i++) {
  160.       yerr[i] = h * h * (y[i] - y0[i]); /* error is third order */
  161.       yim1[i] = y0[i];
  162.     }
  163.  
  164.     /* Make note of step size. */
  165.     state->last_h = h;
  166.  
  167.     return status;
  168.   }
  169. }
  170.  
  171. static int
  172. gear2_reset(void * vstate, size_t dim)
  173. {
  174.   gear2_state_t *state = (gear2_state_t *) vstate;
  175.  
  176.   DBL_ZERO_MEMSET (state->yim1, dim);
  177.   DBL_ZERO_MEMSET (state->k, dim);
  178.   DBL_ZERO_MEMSET (state->y0, dim);
  179.  
  180.   state->primed = 0;
  181.   state->last_h = 0.0;
  182.   return GSL_SUCCESS;
  183. }
  184.  
  185. static unsigned int
  186. gear2_order (void *vstate)
  187. {
  188.   gear2_state_t *state = (gear2_state_t *) vstate;
  189.   state = 0; /* prevent warnings about unused parameters */
  190.   return 3;
  191. }
  192.  
  193. static void
  194. gear2_free(void *vstate)
  195. {
  196.   gear2_state_t *state = (gear2_state_t *) vstate;
  197.  
  198.   free(state->yim1);
  199.   free(state->k);
  200.   free(state->y0);
  201.   gsl_odeiv_step_free(state->primer);
  202.  
  203.   free(state);
  204. }
  205.  
  206. static const gsl_odeiv_step_type gear2_type = { "gear2",    /* name */
  207.   1,                /* can use dydt_in */
  208.   0,                /* gives exact dydt_out */
  209.   &gear2_alloc,
  210.   &gear2_apply,
  211.   &gear2_reset,
  212.   &gear2_order,
  213.   &gear2_free
  214. };
  215.  
  216. const gsl_odeiv_step_type *gsl_odeiv_step_gear2 = &gear2_type;
  217.